#####################################################################
# Table 1 and 2
#####################################################################

rm(list=ls())

library(ggplot2)
library(reshape2)
library(lmtest)
library(sandwich)
library(stargazer)

# load function that does clustered SEs
vcovCluster <- function(
  model,
  cluster
)
{
  require(sandwich)
  require(lmtest)
  if(nrow(model.matrix(model))!=length(cluster)){
    stop("check your data: cluster variable has different N than model")
  }
  M <- length(unique(cluster))
  N <- length(cluster)           
  K <- model$rank   
  if(M<40){
    warning("Fewer than 40 clusters, variances may be unreliable")
  }
  dfc <- (M/(M - 1)) * ((N - 1)/(N - K))
  uj  <- apply(estfun(model), 2, function(x) tapply(x, cluster, sum));
  rcse.cov <- dfc * sandwich(model, meat = crossprod(uj)/N)
  return(rcse.cov)
}

############################
# Read and prepare data 
############################

load("~/Dropbox/Crime Brazil/06_replication/001_rematch_crime.Rdata")

# treatment 
table(d$t_ind)
table(d_match$t_ind)
table(d_rematch$t_ind)

##############################
# Outcome: Iron-fist policies
##############################

# control and fixed effects
mean(d_match$outcome_crime_stronghand_wave4[d_match$t_ind==1])
mean(d_match$outcome_crime_stronghand_wave4[d_match$t_ind==0])
mod1 = lm(d_match$outcome_crime_stronghand_wave4 ~ d_match$t_ind + d_match$crime_stronghand_wave1 + d_match$s6_education_wave2_imp + d_match$s6_education_wave2_miss + d_match$s10_age_wave3_imp + d_match$s10_age_wave3_miss + as.factor(d_match$bairroibge_wave3))
mod1c = coeftest(mod1, vcov = vcovCluster(mod1, cluster = d_match$bairroibge))
round(mod1c,2)

# no control
mean(d_match$outcome_crime_stronghand_wave4[d_match$t_ind==1])
mean(d_match$outcome_crime_stronghand_wave4[d_match$t_ind==0])
mod2 = lm(d_match$outcome_crime_stronghand_wave4 ~ d_match$t_ind + as.factor(d_match$bairroibge_wave3))
mod2c = coeftest(mod2, vcov = vcovCluster(mod2, cluster = d_match$bairroibge))
round(mod2c,3)

# no fe
mean(d_match$outcome_crime_stronghand_wave4[d_match$t_ind==1])
mean(d_match$outcome_crime_stronghand_wave4[d_match$t_ind==0])
mod3 = lm(d_match$outcome_crime_stronghand_wave4 ~ d_match$t_ind + d_match$crime_stronghand_wave1 + d_match$s6_education_wave2_imp + d_match$s6_education_wave2_miss + d_match$s10_age_wave3_imp + d_match$s10_age_wave3_miss)
mod3c = coeftest(mod3, vcov = vcovCluster(mod3, cluster = d_match$bairroibge))
round(mod3c,3)

# no control no fe
mean(d_match$outcome_crime_stronghand_wave4[d_match$t_ind==1])
mean(d_match$outcome_crime_stronghand_wave4[d_match$t_ind==0])
mod4 = lm(d_match$outcome_crime_stronghand_wave4 ~ d_match$t_ind)
mod4c = coeftest(mod4, vcov = vcovCluster(mod4, cluster = d_match$bairroibge))
round(mod4c,3)

###########################################
# Causal mechanism: Support for democracy
###########################################

mean(d_match$outcome_supportdem_wave4[d_match$t_ind==1])
mean(d_match$outcome_supportdem_wave4[d_match$t_ind==0])
mod5 = lm(d_match$outcome_supportdem_wave4 ~ d_match$t_ind + d_match$supportdem_wave1 + d_match$s6_education_wave2_imp + d_match$s6_education_wave2_miss + d_match$s10_age_wave3_imp + d_match$s10_age_wave3_miss +  as.factor(d_match$bairroibge_wave3))
mod5c = coeftest(mod5, vcov = vcovCluster(mod5, cluster = d_match$bairroibge))
round(mod5c,3)

mean(d_match$outcome_supportdem_wave4[d_match$t_ind==1])
mean(d_match$outcome_supportdem_wave4[d_match$t_ind==0])
mod6 = lm(d_match$outcome_supportdem_wave4 ~ d_match$t_ind +  as.factor(d_match$bairroibge_wave3))
mod6c = coeftest(mod6, vcov = vcovCluster(mod6, cluster = d_match$bairroibge))
round(mod6c,3)

mean(d_match$outcome_supportdem_wave4[d_match$t_ind==1])
mean(d_match$outcome_supportdem_wave4[d_match$t_ind==0])
mod7 = lm(d_match$outcome_supportdem_wave4 ~ d_match$t_ind + d_match$supportdem_wave1 + d_match$s6_education_wave2_imp + d_match$s6_education_wave2_miss + d_match$s10_age_wave3_imp + d_match$s10_age_wave3_miss)
mod7c = coeftest(mod7, vcov = vcovCluster(mod7, cluster = d_match$bairroibge))
round(mod7c,3)

mean(d_match$outcome_supportdem_wave4[d_match$t_ind==1])
mean(d_match$outcome_supportdem_wave4[d_match$t_ind==0])
mod8 = lm(d_match$outcome_supportdem_wave4 ~ d_match$t_ind)
mod8c = coeftest(mod8, vcov = vcovCluster(mod8, cluster = d_match$bairroibge))
round(mod8c,3)

##########################
# Tables in LaTex
##########################

# Table 1: Outcome 1
dim(d_match)
stargazer(mod1c,mod2c,mod3c,mod4c,
          type = "text",
          title="Regression results",
          align=TRUE, 
          omit.stat=c("LL","ser","f","rsq","adj.rsq"), 
          no.space=TRUE,  
          multicolumn = TRUE,
          table.placement = "H",
          #single.row = TRUE,
          column.separate = c(2,2),
          covariate.labels=c("Crime victimization"),
          dep.var.caption  = "Strong-handed policies and repression to reduce crime",
          #column.labels = c("Ordinal","Binary"),
          dep.var.labels.include = FALSE,
          #star.cutoffs = c(0.05,0.01,0.001),
          omit = c("crime_stronghand_wave1",
                   "s6_education_wave2_imp",
                   "s6_education_wave2_miss",
                   "s10_age_wave3_imp",
                   "s10_age_wave3_miss",
                   "bairroibge_wave3",
                   "Constant"),
          add.lines = list(c("Controls","Yes","No","Yes","No"),
                           c("Neighborhood fixed effects","Yes","Yes","No","No"),
                           c("Observations","542","542","542","542","542")))

# Table 2: Causal mechanisms
dim(d_match)
stargazer(mod5c,mod6c,mod7c,mod8c,
          type = "text",
          title="Regression results",
          align=TRUE, 
          omit.stat=c("LL","ser","f","rsq","adj.rsq"), 
          no.space=TRUE,  
          multicolumn = TRUE,
          table.placement = "H",
          #single.row = TRUE,
          covariate.labels=c("Crime victimization"),
          dep.var.caption  = "Support for democracy",
          #column.labels = c("Democracy","PT","PMDB","PSDB","PFL"),
          dep.var.labels.include = FALSE,
          #star.cutoffs = c(0.05,0.01,0.001),
          omit = c("supportdem_wave1",
                   "s6_education_wave2_imp",
                   "s6_education_wave2_miss",
                   "s10_age_wave3_imp",
                   "s10_age_wave3_miss",
                   "bairroibge_wave3",
                   "Constant"),
          add.lines = list(c("Controls","Yes","No","Yes","No"),
                           c("Neighborhood fixed effects","Yes","Yes","No","No"),
                           c("Observations","542","542","542","542")))
